Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 20 January 2010 (MN L5T^ style file v2.2) 



Estimates of unresolved point sources contribution to WMAP 5 

L.P.L. Colombo^ & E. Pierpaoli\ 

1 University of Southern California, Los Angeles, CA, 90089-0484 



o 

(N 

o 

(N 



o 

u 

o 



> 

vn 
\o 
rn 

O 

o 



Accepted ???. Received ???; in original form ??? 



ABSTRACT 

We present an alternative estimate of the unresolved point source contribution to the 
WMAP temperature power spectrum based on current knowledge of sources from radio sur- 
veys in the 1.4—90 GHz range. We implement a stochastic extrapolation of radio point sources 
in the NRAO-VLA Sky Survey (NVSS) catalog, from the original 1.4 GHz to the - 100 GHz 
frequency range relevant for CMB experiments. With a bootstrap approach, we generate an 
ensemble of realizations that provides the probability distribution for the flux of each NVSS 
source at the final frequency. The predicted source counts agree with WMAP results for S > 1 
Jy and the corresponding sky maps correlate with WMAP observed maps in Q-, V- and W- 
bands, for sources with flux S > 0.2 Jy. The low-frequency radio surveys found a steeper 
frequency dependence for sources just below the WMAP nominal threshold than the one esti- 
mated by the WMAP team. This feature is present in our simulations and translates into a shift 
of 0.3 — 0.4a in the estimated value of the tilt of the power spectrum of scalar perturbation. 
Us, as well as ujc- This approach demonstrates the use of external point sources datasets for 
CMB data analysis. 
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1 INTRODUCTION 

Modem Cosmic Microwave Background (CMB) observations pro- 
vide a powerful test of our understanding of the Universe. Within 
the generally accepted framework of a ACDM model, the 5- 
year measurements by the Wilkinson Microwave Anisotropy Probe 
(WMAP) constrain the funda mental cosmological parameters with 
a rela tive accuracy of 1 - 10% jPunklev et al.l2009l : lKomatsu et all 
l2009l) . while the upcoming observations by the Planck satellite are 
expe cted to improve these numbers by at least a factor ^ 3 — 4 
(e.g. iThe Planck Collaboratioij l2006l : [Colombo et al. I l2009h . One 
of the main scientific goal s of Planck, and other proposed future 
CMB missions (e.g. EPIC. lBock et"ai]| 20091) , is understanding the 
nature of Inflation. Detection of the B-mode of CMB polariza- 
tion would provide direct evidence of a primordial backg round o f 
Gravitational Waves arising from Inflation ( Kamionkowski et al.l 
Il997l : ISpergel & Zaldarriagalll997h . Even without such detection 
the CMB remains the most powerful probe of Inflation currently ac- 
cessible. A general prediction of inflationary models is that the frac- 
tional amplitude of density fluctuations would be nearly scale inde- 
pendent, so that the corresponding power spectrum could be well 
approximated by an Harrison-Zeldovich form P{k) oc fc"° , where 
the spectral index Us — 1. The amount of deviation ns — 1 con- 
strains the shape of the inflationary potential, and curren t WMAP 
results already allow to rule out several models (iKomatsu et al.l 



l2009h . Significant improvements will come from Planck and the 
next generation CMB missions. 

However, exploitation of the full potential of CMB measure- 
ments requires a deep understanding of instrumental systematics 
and a careful cleaning of foreground contaminants. On small an- 
gular scales, extragalactic point sources are an important source 
of contamination. Bright sources, which can be detected with high 
significance in CMB maps, are typically accounted for by masking 
a small area around the source position during the estimation of 
the CMB angular power spectra, Ce- On the other hand, to a first 
approximation, undetected, and therefore unmasked, point sources 
provide a Poisson noise contribution to the measured Ct as well as 
a non-Gaussiarr signature in the rnaps (e .g. Toffolatti et al.l 



a non-uaussian signature m trie maps (e . g. I lotiolatti et al.l 
'Pierpaoli' '2003'; 'Pierpaoli & Pemj |2004| : ISerra & Cooravl 
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Babich & Pierpaoli 2008). An incorrect determination of this con- 
tamination can lead to relevant biases on the estimated cosmo- 
logical parameters, in particular on fauffenberger et alj|200d 
2008]), whose accurate measurement depends on a careful estimate 
of Ce over the largest range of scales probed by an experiment. In 
this respect the greatest expectations are now on Planck, which is 
a whole sky survey with a small instrumental beam and low noise 
level. If foregrounds and systematics are well under control, Planck 
will be able to confirm or falsify the WMAP findings about Us 
being significantly smaller than one. As inflation predicts that the 
level of departure on ris from one is proportional to the amount of 
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gravitational waves expected, this also has implications on the level 
of B-mode polarization signal expected and is in turn relevant for 
the design and planning of future missions. It is therefore appropri- 
ate to try to approach the issue of point source contamination from 
different angles and conceive different approaches to deal with it. 
In this paper, we undertake this path with an application to WMAP 
data analysis which could lead to new approaches for point sources 
subtraction for Planck. 

At frequencies around ~ 1 GHz there is a wealth of infor- 
mation on extragalactic sources, from nearly full-sky surveys with 
high angular resolution and flux sensitivi ty, like the Green Bank 
Telescope GB6 survey jGregorv et Z|[T99 6). the NRAO-VLA Sky 
Survey (NVSS lCondon et alj|l998h o r the Faint Images of the Ra- 
dio Sky (FIRST iBecker et al.l l2003h . On the other hand, at fre- 
quencies ~ 100 GHz, corresponding to the optimal observational 
window for CMB experiments, there are only few dedicated stud- 
ies covering a handful of sources. Most information comes from 
sources detected in CMB experiments themselves, typically repre- 
senting only the bright end of the point source distribution. In addi- 
tion, this lack of information limited the possibility of extrapolating 
data from low frequencies to the CMB "gold spot" . 

To estimate the residual point sources contribution, the 
WMAP team assumed that unresolved sources followed the same 
frequency scaling as the bright sources detected in the 5 year maps. 
The overall normalization of the residual power spectrum was then 
determined by a fit to maps at three highest frequencies. This ap- 
proach can be considered internal, as it takes into account (mostly) 
only the WMAP data itself, without relying on the low-frequency 
information form other available point sources surveys. As other 
radio surveys provide alternative and further information on ra- 
dio sources properties, incorporating such results in the CMB data 
analysis pipeline could improve the estimation of the extragalac- 
tic sources contamination, or provide an external crosscheck of the 
WMAP results. 

Different authors provided estimates of source counts at CMB 
frequencies on the basi s of the physical properties of the differen t 
source populations (e.g. lToffolatti et al.lll998l : lde Zotti et"al]|2005h . 
However, multifrequency studies, describing the scaling of sources 
from middle frequ encies, v ~ 2Q GH z, to ~ 100 GHz, are now 
startmg to appear jSadleretal.ll2008h and more are expected in 
the future. Together with existing studies describing sources' be- 
haviour from low to intermediate frequencies, these allow to pre- 
dict sources contamination in CMB maps starting from the low fre- 
quency surveys. Such predictions could be used in future data anal- 
ysis to make specific assessments on point sources contamination 
in CMB maps in real space, and account for them with techniques 
different from the Fourier-based ones currently used. 

In this work we follow the empirical approach outlined above 
to assess the unresolved point sources contribution to the WMAP 
data and its impact on cosmological parameters. We use the NVSS 
catalog as a template of source fluxes and positions at radio fre- 
quency, and extrapolate that information to CMB frequencies. At 
difference with previous works, we extrapolate each individual 
source in the NVSS rather than generating source populations at 
CMB frequencies by randomly sampling a theoretical distribution, 
or by extrapolating only the statistical properties of the source pop- 
ulation (e.g. Waldram et al. 2007). This work therefore provides a 
crosscheck of WMAP estimates of undetected point sources con- 
tamination, and hence of the determination of and other param- 
eters. As mentioned above, since this approach maintains the infor- 
mation on the actual spatial distribution of sources, it can be useful 
in cleaning future Planck maps or to estimate the impact of cluster- 



ing. Here, however, we refrain from using an alternative likelihood 
approach to the treatment of point sources in the parameter estima- 
tion and adopt the WMAP strategy to that aim. 

The outline of the paper is as follow. In Section |2] we briefly 
review the contribution of unresolved sources to CMB temperature 
spectra. Section [3] describes our extrapolation procedure of NVSS 
source to WMAP frequencies. In section |4] and |5] we discuss the 
implication of our approach at the number counts and maps level, 
respectively, while in section|6]we debate the corresponding impact 
on determination of cosmological parameters. Finally, in section |7] 
we draw our conclusions. 



2 UNRESOLVED SOURCES CONTRIBUTION TO CMB 
SPECTRA 

The angular resolutions of WMAP (> 13 arcmin) and Planck 
(> 5 arcmin) are significantly larger than the angular dimen- 
sions of typical extragalactic radio sources, which can then be ef- 
fectively considered as point-like objects in the resulting CMB 
maps. Several works studied the secondary contribution to CMB 
anisotropi es due to point sources below the detection t hresh- 
old (e.g. If ranceschini et al. 1989; Tegm ark & Efstathioul 1 1 9961 : 
iToffolatti et al. 1998; Scott & White 199^. In this section we 
briefly summarize the main results. 

Unresolved sources contribution to Ci can be divided into a 
Poisson term due to the random sources positioning, and a correc- 
tion accounting for clustering in the sources distribution. We con- 
sider CMB observations at a single frequency and suppose that 
all sources above a limit flux Sc will be detected and masked from 
the final map, and assume a source population characterized by the 
differential counts dN{> S)/dS. 

The Poisson noise contribution to the measured CMB Ci is 
then given by: 



Ci =g {vo) 
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where g{i') converts from flux density to thermodynamic tempera- 
ture: 
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The source correction does not depend on £ and, since approx- 
imately Ci oc at high multipoles up to £ ~ 2000, the 
source correction becomes relevant on angular scales smaller than 
180/^ ~ 180/1000 ~ 0.2 deg. The correction due to clustering is 
given by: 
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here uje is the harmonic transform of the source two-point angular 
correlation function. 

Real CMB experiments typically have complicated noise 
properties, due to a combination of effects including instrumental 
systematics, scanning strategy and foreground removal. As a conse- 
quence, the probability of detecting a point source with flux So will 
not be uniform on the whole sky and Sc will not be well defined. 
This is not irrelevant as the unresolved point sources spectrum in 
eq[T]is typically dominated by the strongest amongst the residual 
sources. While the source population above > 1 Jy in the WMAP 
5 year catalog ( WM APS, i Wright et al. 200 9) is well described by 
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Figure 1. Fraction of unresolved source contribution to WMAP spectra 
attributed to sources with fluxes St., assuming that all sources with 
S ^ 0.7 Jy have been removed from the maps. We used the de Zotti et 
al.(2005) model for the number counts of sources with S < 0.7 Jy. Sources 
below 0.1 Jy contribute < 10% of the total unresolved power. 

an Euclidean scaling of the number counts dN/dS oc S~'^'^ , de- 
tected fainter sources suffers from significant completeness issues. 
Depending on the frequency band considered, > 90% of detected 
sources have a flux ^ 0.7 Jy. Ass uming then an optim istic detec- 
tion threshold = 0.7 Jy and the lde Zotti et al.l l l2005h model for 
dN/dS, we can estimate that sources with 5* < 0.1 Jy contribute 
< 10% of the total unresolved source power in WMAP data, as 
shown in figure [T] On the contrary, for Planck, whose 90% detec- 
tion confidence lim it is expected to be Sc ~ 0.2 Jy at 100 GHz 
jVielva et alj|2003h . sources with S ~ 0.1 Jy will play a dominant 
role in determination of C^^". Therefore, the accuracy required for 
the extrapolations at a given flux level depends on the characteris- 
tics of the target experiment. 



3 EXTRAPOLATION PROCEDURE 
3.1 Reference Catalogs 

Instead of extrapolating sources directly from NVSS to WMAP 
frequencies, we consider a number of intermediate steps, depend- 
ing on the initial flux of the source considered. This is particularly 
relevant at frequencies below 20 GHz, where a singl e power law 
is not an accurate description for many sources (e.g. I Sadler et al.l 
[2006). By introducing several intermediate steps, we account for 
the diversity of behaviours observed in the data. The number and 
frequency range of the steps used here was determined by avail- 
able data. In selecting the reference catalogs that we will use for 
the extrapolation, our goal was to get a coverage of fluxes relevant 
for the estimation of unresolved sources contribution to WMAP 
data over the full range of frequencies considered. While this is 
not an issue in the ~ GHz range, coverage becomes sketchier and 
patchier at higher i/. According to the discussion in section [2] the 
main contamination to WMAP spectra is due to undetected sources 
with fluxes 0.1 ~ 0.7 Jy in the CMB frequency range (40 — 90 



GHz). We will show this population to be dominated by sources 
with flat spectrum (q > —0.5 and 5* oc u"), therefore correspond- 
ing to sources with fluxes from a few hundreds of mjy to a few Jy 
at 1.4 GHz. 

In order to perform the extrapolation, we considered the fol- 
lowing sets of measurements: 

• The NRAO-VLA Sky Survey (NVSS): it covers the full sky 
north of —40 deg declination with an angular resolution of 45" 
at 1.4 GHz. It includes about 1.8 x 10^ discrete sources with a 
completeness limit of 5*1.4 = 2.5 mJy. 

• The Green Bank Telescope Survey at 6cm (GB6) observed the 
sky in the declination band deg < S < 75 deg with an angular 
resolution ~ 3.5' (however, the actual beam is not circular). It in- 
cludes 75162 sources above 5*4.8 = ISmJy. 

• The NVSS followup measurements bv lMason et al.l ( |2009|) . A 
set of 3165 NVSS faint sources with 5 i.4 < 100 mJy and f alling 
in the Cosmic Background Imager (CBl lReadhead et al.l20Q4h field 
was observed atv = 31 GHz using the Green Bank Telescope and 
the Owens Valley Radio Observatory. The resulting set of observa- 
tions allowed to derive the distribution of spectral indexes of for 
sources with 5i.4 ^5 0.1 Jy. 

• The Ninth Cambridg e Surv ey at 15 GHz (9C) followup obser- 
vations bv fWaldram et alj ( l2007h . The catalog includes 121 sources 
selected at 15.2 GHz with a completeness limit of 5i5 = 25 mJy 
and simultaneous measurements at 4.85, 15.2, 22 and 43 GHz. As 
the 9C survey covers the sky area observed by the Very Small Ar- 
ray (VSA Watson et al. 2003), which was selected to be devoid of 
bright sources, the sample contains only a handful of sources with 
5 > 100 mJy. 

• The AT20G followup observations bv lSadler et allfeOOSl) . The 
work presents two samples of sources selected at 20 GHz with si- 
multaneous measurements at 90 GHz: a first set of 59 inverted spec- 
tra sources with S20 > 50 mJy, and a flux-limited sample compris- 
ing 70 sources with 52o > 150 mJy. We consider here only the 
flux-limited sample, which constitutes our main reference sample 
for the extrapolation of sources with 52o ^ 1 Jy from 20 to 94 
GHz. 

• At frequencies above ~ 23 GHz and for fluxes above ~ 1 Jy, 
the most comprehensive survey of point sources is provided by the 
WMAP 5 year catalog. WMAP observed the sky at five broad bands 
K, Ka, Q, V and W with central frequencies of ~23, 33, 41, 61, 94 
GHz respectively. The catalog includes a total of 390 sources which 
have been detected with a 5a confidence level in at least one of the 
bands, and is complete for fluxes above ~ 1 Jy in all bands. 

3.2 Method 

Since the number of sources with high frequency information is 
smaller than the number of NVSS sources, reconstructing the ac- 
tual scaling of individual sources is not possible. Our assumption is 
thus that from the reference catalogs discussed above, we can build 
reference samples which are representative of the behaviour of the 
whole source population in the respective range of frequencies and 
fluxes. We then perform a set of bootstrap simulations in which at 
each extrapolation step, a simulated source is randomly paired to a 
reference source in the catalog covering the next higher frequency 
range considered, and we use the measured spectral index of the 
reference source to further extrapolate the simulated source. Some 
of the catalogs considered provide multifrequency information for 
their sources. When this is the case, we use the set of spectral in- 
dexes of the catalog source to extrapolate fluxes over the range of 
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frequencies considered. For each source in the NVSS catalog, we 
generate a set of 800 simulations. While in two different simula- 
tions a given source can be extrapolated in different ways, the set of 
simulations provides the probability distribution for flux at WMAP 
frequencies. 

In broad terms, we define two main frequency ranges: 1) the 
interval from NVSS to the lowest WMAP band, 1.4 - 23 GHz, 
and 2) the WMAP frequency range, 23 — 94 GHz. Sources with 
'S' i.4 100 mjy are extrapolated directly to 23 GHz according 
to lMason etal] j2009l) . while for those with Si a ^ 300 mJy we 
intr oduce a number of int ermediate steps based based on the GB6 
and lWaldram et al.l ( l2007h measurements. In the intermediate flux 
range, we weight the two approaches in order to fit the 33 GHz 
low flux counts. From 23 GHz upwards, sources with 523 > 1 Jy 
are propagated based on the WMAP5 catalog, while for lower flux 
sources we base the extrapolation on the Sadler etal^ C2008.') and 
IWaldrametZli2007h reference samples. 

In the following, we discuss each step in detail: 

• The 1.4 GHz to 4.85 GHz reference sample. 'Kimba ll & Ivezid 
( I2OO81) provide a unified catalog of sources from four radio surveys 
NVSS, GB6, FIRST, WENS and from the optical SDSS surveys, 
matching sources in the different catalogs. From the NVSS-GB6 
pairings identified in that work, we construct a reference sample 
of sources which provides spectral indexes between 1.4 GHz and 
4.85 GHz. We consider only pairs in which the GB6 counterpart 
is found within 70" of the NVSS source, corresponding to an es- 
timated completeness of 0.979 and an efficiency (fraction of de- 
tected pairs corresponding to actual physical matches) of 0.79 (for 
further details see lKimball & Ivezidl2 008). In addition, we exclude 
sources flagged as extended in the GB6 survey and include only 
unique pair of sources, i.e. sources in one catalog with multiple 
matches in the other catalog are excluded. Since multiple sidelobes 
of a single source may show up as different NVSS components, in 
this way we may bias the pairs toward compact sources. However 
NVSS resolution is such t hat less than 1% of sources is resolved 
into multiple components jBest et al]|2005h . compared with ~ 7% 
of GB6 sources with multiple NVSS matches in the reference cat- 
alog. We then conclude that the majority of multiple matches we 
find are spurious pairs. 

The distribution of spectral indexes so obtained describes the 
1.4—4.85 GHz behaviour of sources selected at the GB6 frequency, 
and can be used to accurately propagate sources from 4.85 GHz to 
the lower frequency. In order to use this distribution to extrapo- 
late from 1.4 GHz to higher frequencies, some adjustments are re- 
quired. The NVSS survey has a flux limit of 2.5mJy, while the GB6 
includes sources above 18mJy, therefore low-flux NVSS sources 
which appear also in the GB6 survey will be biased to have a ris- 
ing spectral index. Figure |2] shows the average a* ;| as a function 
of Si. 4,. While for S-1.4, > 0.3 Jy the spectral index is a decreasing 
function of the NVSS flux, for 51.4 < 0.1 Jy af;! shows a steep 
upturn due to completeness issues. We therefore exclude from the 
reference sample all sources with 5'i.4 ^0.1 Jy. We further divide 
the reference sample in 10 logarithmic bins covering the flux range 
0.1 — 3.0 Jy, depending on the value of 51.4, with an additional 
bin including all sources with S-1.4, > 3 Jy. In the end the reference 
sample includes 31000 sources, with corresponding spectral index 
Qi and uncertainty Soti. 

In principle, since the reference sample is based on sources with 
matches in both NVSS and GB6 catalogs, we may be missing: a) 
sources with strongly rising spectra, which are absent in NVSS but 
are detected by GB6, or b) sources with very steep spectra, which 




s,., (Jy) 



Figure 2. The average 1.4 — 4.85 GHz spectral index of GB6 sources as 
a function of 1.4 GHz flux. Below 5i,4 ~ lOOmJy NVSS sources with a 
GB6 counterpart are biased to have upturn spectra. 

fall below the GB6 detection threshold. IWaldrametid] ( l2009l) find 
that 4.3% of sources selected at 15.2 GHz with lOmJy ^ ^is ^ 
15mJy do not have an NVSS counterpart. GB6 sources are se- 
lected with hig h er thre shold (ISmJy) and at lower frequencies than 
IWaldram et aU | |2009|) . Therefore, sources of type (a) would re- 
quire spectral indexes significantly more rising than those of the 
rising population discussed in that work, an we expect the frac- 
tion of sources detected at 4.85 GHz but not in the NVSS to be 
significantly lower than IWaldram et all ( |2009|) estimates. Regard- 
ing point (b), we note that sources in the reference sample with 
100 5*1.4 ^ 200 mJy have an average spectral index af;! — 
—0.84 ±0.34. Therefore we estimate the fraction of NVSS sources 
with no GB6 counterpart to be ~ 5% at lOOmJy and ~ 0.5% at 
150 mJy. 

• From 1.4 GHz to 23 GHz for sources with Si. 4 ^ 300 mJy. 
From the previous discussion, we conclude that the thinned ref- 
erence sample defined above provide an accurate description of 
the frequency scaling from 1.4 to 4.85 GHz for the population of 
sources with S1.4 > 300 mJy, which are responsible for the bulk of 
source contamination at CMB frequencies. Therefore, each NVSS 
source with 5*1.4 > 300 mJy is randomly paired to a spectral in- 
dex Qi from the corresponding flux bin of the thinned reference 
sample. We then randomly draw a spectral index a from a Gaus- 
sian with mean a; and width Sai. In this way we incorporate in 
the extrapolation the measured uncertainty on ai, however, since 
the distribution of the reference spectral indexes is not necessar- 
ily Gaussian, the overall distribution of the a both within a single 
realization and among different realizations will not be Gaussian. 
The flux of the NVSS source is extrapolated to 4.85 GHz assuming 
S ct 1/°' . The source is then propagated to 23 GHz using the spec- 
tral indexes by Waldram et al. (2007). Since the flux threshold of 
the 15.2 GHz sample is higher than the GB6 limit, we need to thin 
the catalog to avoid that sources with low fluxes at 4.85 will bias 
the sample toward flatter spectral indexes. Due to the low number of 
sources, we cannot study the behaviour of ct{S) as done in the pre- 
vious section. Instead we compute the average spectral index 
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using all the 121 sources in the sample. Using this average a\%, 
we extrapolate the 15.2 GHz threshold to 4.85 GHz and remove 
form the 4.8 — 15.2 GHz reference sample all sources with 5*4.8 
below this value. The process is iterated until no more sources are 
rejected. We are left with 76 sources with S'4.85 > 75 mjy, which 
we assume provide a representative description of the behaviour of 
sources between 4.85 and 23 GHz. Each source simulated at the 
previous step is paired with a random source in this new reference 
sample, and its corresponding set of spectral index. We use this set 
to propagate the source to 23 GHz, going through an intermediate 
step at 15.2 GHz. 

• From 1.4 GHz to 23 GHz for sources with 6*1.4 ^ 100 mJy. 
For low flux sources, 51.4 < 100 mJy, we extrapolate directly to 
23 GHz assuming S oc v", wi th a randomly draw n from the dis- 
tribution of spectral indexes bv lMason et al.l j2009l) . In order to be 
consistent with the treatment of sources with Si. 4 > 100 mJy de- 
scribed above, we assumed that such distribution be valid in the 
range 1.4 — 23 GHz, even if it is based on 31 GHz foUowup of 
NVSS sources. 

• From 1.4 GHz to 23 GHz for sources with lOOmJy < 5*1.4 < 
300 mJy. Sources with lO OmJy < 5i.4 < 300 mJy are above the 
upper validity limit of the iMason et al.l j2009l) study, and may yet 
be slightly affected by spurious effects in the NVSS - GB6 refer- 
ence sample defined above. Therefore, in this intermediate regime, 
we randomly chose between the two approaches discussed above, 
and tune the selection function to reproduce low flux counts at 
higher frequencies. The probabi lity of extrapolating sources using 
the spectral index distribution of iMason et al.l | |2009|) is given by: 

5i.4 — St 



P(5l.4 



tanh 



(4) 



The values 5, and A are chosen by fitting the 33 GHz num- 
ber counts from the extrapolation procedure described here to the 
observed number counts by DASI and VSA in the 0.05 — 1 Jy 
range. Sources with lower fluxes are not a significant contaminant 
in WMAP maps, as discussed above. Note that we do not require 
the selection function in equation|4]to be continuous at the lOOmJy 
and 300mJy boundaries. 

• From 23 GHz to 94 GHz for sources with S23 ^ 1 Jy. Sources 
with 523 > 1 Jy are extrapolated to higher frequencies accord- 
ing to the WMAP5 catalog itself. We divide the WMAP5 cata- 
log in two sub catalogs comprising sources with 823 > 2 Jy and 
IJy < 523 < 2 Jy, respectively containing 69 and 166 elements. 
When a simulated source has S23 > 1 Jy, we randomly pair it with 
a spectral index from the subcatalog covering the corresponding 
flux range, and extrapolate the source to higher frequency using a 
simple power law. 

• From 23 GHz to 94 GHz for sources with S23 < 1 Jy. 
For fluxes below ~ 1 Jy the WMAP catalog is not complete. 
IWaldrametall feOO?) suggest a procedure to extrapolate fluxes into 
the WMAP frequency range based on the spectral behaviour of the 
source around ~ 20 GHz. For sources for which all ^ f^ii ^5 
they suggest fitting a quadratic form to the log(5) — log(z/) relation, 
while for the remaining sources they adopt a linear log (5) — log 
using 022- Alternatively, we co nsider the set of spectral indexes 
based on the lSadler et al.l j2008i) sample. Taking into account both 
sets of observations, we consider then the following extrapolation 
strategy for sources with S23 < 1 Jy: 

- Sources with 523 < 100 mJy are randomly paired to source 
of IWaldram et al.l ( l2007h and extrapolated according to the pro- 
cedure described there; 

- Sources with 523 > 400mJy are propagated assuming a 



single power law with spectra index d rawn randomly from the 
reference sample of I Sadler et alj ( l2008h : 

~ Sources with 100 < S23 < 400mJy are extrapolated 
choosing randomly between the previous two methods; the prob- 
ability of selecting one method over the other is given by the 
relative number of sources in the relevant flux range. 

In addition, when at a given step we have different extrapo- 
lation procedures depending of the flux of the source, we do not 
switch abruptly from one regime to the other, but we adopt a lin- 
ear transition between the two with a width which is the greater of 
lOOmJy and 10% of the threshold flux. Moreover, we do not con- 
sider polarization. 



4 RESULTS 

4.1 Flux distribution 

This approach provides both the probability that a source with an 
NVSS flux 5i.4 will have a 94 GHz flux S94, i.e. P(594|Si.4), as 
well as the complementary distribution P(5i.4j5g4). In figure[3]we 
plot the distribution of the initial 1.4 GHz fluxes for sources with 
final 94 GHz fluxes of S94 ~ 0.01, ~ 0.10 and ~ 1.0 Jy, which 
provides an estimate of P(5i.4|594). The plots refer to the full set 
of 800 simulations. Figure|4]shows instead P(594|5i.4). A general 
expectation is that for flux levels relevant to CMB experiments, the 
high frequency population of sources is dominated by objects with 
flat spectral indexes. We recover this behaviour in the simulations, 
as shown in figures [3] and |4] In particular, according to the discus- 
sion of section |2l the residual contamination in WMAP maps will 
be dominated by sources with 5 > 0.1 Jy in the frequency range 
60 — 90 GHz, and we expect the majority of these sources to have 



5i 



lOOmJy. Inaccuracy in the extrapolations of sources with 



lower 5i.4 will have a minimal impact on estimates of the residual 
contamination in WMAP maps. 



4.2 Number counts 

For each simulation, we compute the corresponding differential 
number counts dn{S)/dS and then average over the whole set of 
simulations. In figure [5] we compare the incomplete counts from 
WMAP5 source catalog with the ensemble average differential 
counts at the frequencies of 33, 41, 61 and 94 GHz, approximately 
corresponding to the central frequencies of WMAP Ka, Q, V and 
W bands. At 33 GHz we compare our estimation also with results 
from CBI, VSA, and DASI sur veys, while at 94 GHz we also plot 
predi ction from earlier works jde Zotti et alJl200^ ; IWaldram et al.l 
I2OO7I) . 

Comparison with lower flux data at 33 GHz, shows that the 
methods discussed here tends to underestimate source counts be- 
low 5 ~ 0.02 — 0.03 Jy. These sources do not provide a relevant 
contribution to the total residual contamination to WMAP spectra, 
but will play a more relevant role for Planck data analysis. Since the 
uncertainty in our methods is dominated mostly by the low number 
of reference sources above 4.85 GHz, we expect a better agreement 
as the pool of reference sources increases. In the range 0.5—1.0 Jy, 
simulations show evident artifacts due to the transition between dif- 
ferent regimes of flux extrapolation, i.e. using spectral index from 
the WMAP5 catalog for S23 > 1 Jy. They are due to the low num- 
ber of sources in that flux range in the base catalogs; as above, we 
expect these artifact to significantly decrease as more data over the 
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Figure 3. Distiibutions of the initial 1.4 GHz fluxes for sources with a final 94 GHz as indicated in the respective panels. Histograms refer to the full set of 
800 simulations and have been normalized to have unit area. 
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Figure 4. As in figure[3]but showing distributions of simulated 94 GHz fluxes for sources with initial 1 .4 GHz flux as indicated in the different panels. 



relevant frequency and flux range will be included in the reference 
catalogs. 

While our extrapolations correctly recover WMAP5 results at 
lower frequencies, counts at 94 GHz are overestimated by ~ 50%. 
In figure[6]we plot the integrated counts A^(> S — IJy) from our 
simulations and from WMAP5 data. For reference, we also show 
an analytic e stimate of the coun ts assuming dN/dS oc S~^'^ and 
5* oc u'^ ^^ dWright et alj2009h . normalized to fit WMAP5 counts 
in K, Ka and Q bands. WMAP5 data show a steepening of the 
counts in W band compared to lower frequencies, while our counts 
seem to be in better agreement with a linear log(5') — log(!^) re- 
lation. As discussed in section [X2l our extrapolations of sources 
with S23 > 1 Jy is mainly based on the family of spectral in- 
dexes provided by the WMAP5 team, which can be approximated 
by a Gaussian distribution with mean a = —0.09 and dispersion 
(Tq = 0.17. Therefore, as expected source counts in the simula- 
tions are in agreement with the analytic estimate plotted in figure 
[SJ the slightly steeper behaviour seen in simulations is compati- 
ble with the actual distribution of WMAP5 spectral indexes hav- 
ing a tail toward steeper a. This suggests that there is some ten- 
sion between the WMAP5 W-band counts and their spectral index 
distribution. A possible explanation could be a progressive steep- 
ening of the spectral index of bright sources with increasing fre- 
quency. [cionzalez-Nuevo et al. (2008) noticed that the spectral in- 
dex of WMAP sources in the 41 — 64 GHz range is steeper than 
their spectral index in the 5 — 23GHz range, and radio source mea- 
surements in the ~ 100 — 250 GHz range by the QUaD tele- 

" ! 

Friedman et al. 2009) and the South Pole Telescope (SPT, 
IVieira et al..2009.) provide additional evidence in this direction (see 



Ide Zotti et a l.'2009*. for a recent review). If this is actually the case, 
describing the source behaviour in the 23 — 94 GHz range with a 
single power law may not be entirely accurate. 

Alternatively, there may be some unaccounted for systematics 
affecting WMAP detection efficiency in W band. The WMAP team 
required that for a source to be included in the catalog it needed to 
be detected at more than 5a, in at least one band. Its flux in the other 
bands would be included if: a) it was measured at more than 2a, b) 
the fitted source w i dth is within a factor of 2 of the actual beam 
jWright et al]|2009h . lSawangwit & Shaiiki feOOa") pointed out that 
the W band beam response to discrete sources appears to depend in 
a non-linear fashion on source flux. In particular, they find that the 
effective beam profile is significantly wider for point sources than 
for the Jupiter observations used as a basis for the WMAP team's 
beam analysis, which can create a significant bias in determining 
if a source fits requirement (b). Since number counts from simu- 
lations do not include systematics due to WMAP actual detection 
procedure, this effect, if confirmed, could be an alternative expla- 
nation for the discrepancies between simulated and actual number 
counts. Further data on source behaviour at high frequency, e.g. 
from Planck, will help solve this issue. 

In addition, all WMAP sources with measured W band flux 
have been detected also in at least one of the lower frequencies. 
As long as a WMAP source is detecte d in at least one b and, it is 
masked in the power spectrum analysis dNolta etakll 20091) . and we 
follow the same procedure in building sky mask for the analysis 
of section |6l Therefore, any mismatch in the number of W band 
bright sources between simulations and real data will not impact 
our estimates of the unresolved source power as long as counts at 
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Figure 5. Compaiison of the predicted Euclidean counts from tlie ensemble of NVSS simulations of this work (solid lines) with the incomplete counts from 
the WMAP 5 year catalog (squares), at the frequencies of 33, 41, 61 and 94 GHz. At 33 GHz we also plot data from the CBI, VSA and DASI survey, while at 
94 GHz we compare with predictions from earlier works. 



the lower frequency for the simulations are in good agreement with 
WMAP findings. 



5 MAPS 

An advantage of this approach consists in that it naturally provides 
an estimates of bright sources position and accounts for cluster- 
ing of unresolved sources. However, agreement of number counts 
estimated from simulations with those computed by data does not 
guarantee that the method correctly recovers the spatial distribu- 
tion of sources. In order to check whether the simulations match the 
structure observed in WMAP maps, for each realization we produce 
source-only (i.e. without noise and CMB) maps and compare them 
with WMAP observed maps. Source maps also allow and inde- 
pendent determination of unresolved source contamination. Since 
WMAP estimate of this contribution is based on the Q, V and W 
bands data, we generate map for each of the ab ove frequencies. 
We place the sources on an HEALPIjfl dGorski et al.1l2005l) 



A 



h I 



WMAP5 
Simulations 
3 xy-"' 



K 



Ka 



V 



W 



Figure 6. Comparison of integrated counts in simulations (red/triangles) 
and in WMAP 5 year catalog (black/squares). The dashed line shows the 
expected behaviour assuming a scaling of differential counts dn/dS oc 
S^^'^' and a frequency dependence 5 oc u—0 09 jj^g curve amplitude is 
fixed by a fit to WMAP count in K, Ka and Q bands. In computing the 
number counts we excluded the sky area with declination ^ — 40 dog and 
within the WMAP Syr point source catalog mask. 



http://healpix.jpl.nasa.gov/ 
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map at a resolution Nsidc = 2048, corresponding to a pixel size of 
~ 1.7'. This value is considerably greater than the NVSS beam and 
angular resolution, so we can consider the sources as points when 
placing them on the map, and considerably smaller than WMAP W 
beam of ~ 13', so that any artifacts due to pixelization, smoothing 
etc. will be minimal. The actual WMAP beams are not Gaussian, 
their exact profile varying between each set of Differential Arrays 
(DAs). To reduce the number of maps we have to deal with, we 
consider individual frequencies, rather than individual DAs. Each 
map is then smoothed with an effective beam given by the average, 
in harmonic space, of the beam profiles of the DAs of the corre- 
sponding frequency. The smoothed maps are then degraded to a 
resolution Nsidc = 512. 

In order to check that our predictions recover the point source 
structure seen in the maps, we compute the linear correlation co- 
efficient between each individual simulation and the corresponding 
WMAP Syr map: 

^ ^ ((X, - (X)) (M, - (M)))5>5, ^ 

Here X and M respectively refer to the simulated and WMAP 
maps, after application of an harmonic space filter bi/{b'^Ce + 
(7"°"^"), where be is the harmonic transform of the beam profile, 
Ce is a th eory CMB power spe ctrum based on WMAP5 best fit pa- 
rameters tounklev et al. 200?) and 01°^^" is the noise power spec- 
trum j Wrig ht et al. 2009). The means and variances of the maps, 
(X) and (Tx, and the correlation coefficients are computed using 
only pixels within radius equal to the nominal frequency FWHM, 
of sources with flux above a threshold Sc- To assess the signifi- 
cance of this result, we compare the value of r computed from the 
simulations-WMAP5 correlations, with that computed by correlat- 
ing our simulations with a mock microwave sky which includes 
contribution from CMB, WMAP5-like noise, Galactic foregrounds 
and a point sources population obtained according to the proce- 
dure described in this work. In practice, we randomly chose one 
of the simulations as the actual distribution of point sources, and 
add to this a CMB realization, an anisotropic Gaussian noise term 
and a foreground contribution. The pixel noise variance is given 
by n = no/\/7VobI, where A'^obs is the effective number of ob- 
servations in each pixel and np = 2 . 197 mK for the Q band and 
6.547 mK in W band jLimonet al. 1120091) . However, we do not 
account for noise correlations between different pixels. The fore- 
ground contribution is based on the Maximum Entropy Method 
maps0 (Gold et al. 2009). The maps are provided at a resolution 
Nsidc = 128 and smoothed with a 1 deg beam, and are thus have 
no power at scales £ > 300. Therefore the foregrounds contribution 
to the variance and mean of the mock sky, in particular to the small 
scales relevant to point sources, will be significantly lower than the 
impact on the actual WMAP sky maps. We generate different sets 
of mock skies, corresponding to different point sources, CMB and 
noise realizations. The mock sky then constitute a best case sce- 
nario, in which the point source distribution follows exactly the 
model discussed in this work, and several systematic effects, like 
foregrounds and correlated noise, are absent or play a reduced role. 
Tables [T] and |2] show results of the comparison in Q and W bands, 
respectively. 

A general remark is that even in the best case scenario we 
expect only a moderate level of correlation: r < 0.40 and r < 0.25 
for bright sources, respectively in Q and W band. As expected, the 



^ http://lambda.gsfc.nasa.gov/product/map/dr3/mem_mapsjnfo.cfm 



degree of correlation between simulations and actual maps is lower 
than that between simulations and mock data. While for W band the 
difference is within the associated standard deviation, in Q band the 
discrepancy is of order 2a for bright sources increasing to ~ 3cr for 
faint ones. Foreground contamination is more relevant at 41 GHz 
than at 94 GHz and this may account for the significant loss of 
correlation in Q band. 



6 POWER SPECTRA AND EFFECT ON PARAMETERS 

The WMAP team estimated the amplitude of the unresolved point 
source contribution to a cross spectrum Cj^^^ (Y1Y2 represents 
a pair of different DAs, e.g. Q1Q2, QiVi) assuming a power-law 
scaling in frequency: 

cr = Ao,K).(^.)(^)"(^)", (6) 

where ui, 1/2 are the frequencies of the DAs considered, g{i') was 
defined in equation |2l and the spectral index a is assumed to be 
either 0, as most sources are expected to have a flat spectrum, or 
-0.09, corresponding to the mean of of a Gaussian fit to the spec- 
tral indexes of the actual WMAP5 catalog I Wright et alj2009t ). The 
common amplitude Aq = 0.011 ± 0.001 fj,K^sr is determined by 
fitting the spectral shape of equation |6] to the cross spectra from 
the Q,V,W bands; the estimated Ap has only a sma ll dependence on 
the two values of a considered jNolta et alhoogi) . The final power 
spectrum, instead, is a linear combination of the possible VV,VW, 
and WW cross sp ectra, weighted by the corresponding inverse 
covariance matrix l lHinshaw etU] l2003l : iNolta et al.l 120091) . This 
method is internal, as it includes information only from WMAP 
own measurements. On the contrary, the approach discussed in 
this paper, provides a partially external method, in the sense that 
it incorporates both information from WMAP data (to propagate 
sources with 5*23 > 1 Jy) and from other surveys, and could pro- 
vide a cross check of WMAP results. 

To estimate the unresolved point sources contamination from 
our simulations we need to define an individual mask for each 
realization. We start by applying the same Galactic cut as for 
the KQ8 jfl mask, used in the WMAP5 power spectrum analysis 
jNolta et al]|2009l) , and also mask the area outside of the NVSS 
sky coverage. This base mask is the same for all realizations. In 
addition, the KQ85 mask removes the sky area within 0.6 deg of 
sources in the WMAP5 catalog. The WMAP detection efficiency 
for radio sources depends on a number of factors, including fre- 
quency, position on the sky and the actual CMB and noise fluc- 
tuations around the source. Replicating all these effects in all our 
simulations would be too time consuming and instead we adopt a 
simplified procedure to "detect" and mask sources. The WMAP5 
catalog is complete for sources > 1 Jy at all frequencies, and the 
number counts for those sources follow an Euclidean distribution. 
Therefore, at each frequency we assume a power law shape for the 
number counts F^S~^'^ , and fix the value of by fitting to the 
observed WMAP5 differential number counts nv{S) for S > 1 
Jy at the appropriate frequency. We then assign to each simulated 
source a detection probability P{S,) = n^{S:,) / {F^S-^''). On 
average, in each simulation we detect 304 ± 16 sources, compared 
with 321 actual WMAP5 sources in the same sky area, and we mask 
a 0.6 deg radius around each detected source. Alternatively, we 

^ http://Iambda.gsfc.nasa.gov/product/map/dr3/masks_info.cfm 
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Sc (Jy) 


WMAP5 


mock 


mock 


mock 


mock 


mock 


mock 




maps 


data 1 


data 2 


data 3 


data 4 


data 5 


data 6 


1.0 


0.22 ±0.07 


0.34 ± 0.08 


0.36 ±0.07 


0.37 ±0.08 


0.35 ± 0.07 


0.32 ±0.07 


0.37 ±0.07 


0.5 


0.18 ±0.05 


0.30 ±0.06 


0.32 ±0.05 


0.31 ±0.06 


0.29 ± 0.05 


0.26 ±0.05 


0.31 ±0.05 


0.2 


0.13 ±0.04 


0.23 ±0.04 


0.25 ± 0.04 


0.23 ±0.04 


0.22 ± 0.04 


0.19 ±0.03 


0.25 ± 0.04 


0.0 


0.03 ±0.02 


0.06 ±0.02 


0.07 ±0.01 


0.06 ±0.01 


0.06 ± 0.01 


0.05 ±0.01 


0.08 ±0.01 



Table 1. Linear con'elation coefficient r between simulated source maps at 41 GHz and WMAP 5 year Q band maps, for different flux cuts Sc- Only pixels 
within IFWHM radius of sources with 5 > 5c are included in the computation of r. For comparisons, columns 3-8 show the same quantity computed by 
correlating our set of simulations with mock WMAP-like skies for which the point source component follows the same distribution as our simulations. 



Sc (Jy) 


WMAP5 


mock 


mock 


mock 


mock 


mock 


mock 




maps 


data 1 


data 2 


data 3 


data 4 


data 5 


data 6 


1.0 


0.14 ±0.08 


0.21 ±0.09 


0.24 ±0.08 


0.20 ±0.08 


0.18 ±0.07 


0.16 ± 0.07 


0.18 ±0.07 


0.5 


0.14 ±0.06 


0.18 ±0.07 


0.21 ± 0.06 


0.16 ±0.06 


0.15 ±0.05 


0.14 ±0.05 


0.15 ±0.05 


0.2 


0.12 ±0.04 


0.13 ±0.04 


0.16 ± 0.04 


0.11 ±0.03 


0.11 ±0.03 


0.11 ±0.03 


0.11 ±0.03 


0.0 


0.014 ±0.003 


0.011 ±0.003 


0.014 ±0.003 


0.010 ±0.002 


0.010 ±0.002 


0.010 ±0.003 


0.011 ±0.002 



Table 2. As tableHbut for W band. 



consider the conservative approach of masking only sources with 
a Q band flux S41 > 1.0 Jy. 

We apply the resulting mask to the corresponding set of maps 
and compute t he QQ, VV, VW and WW spectra using a MASTER 
approach (Hiv on et alj|2002l). The VV. VW and WW spectra are 
then combined according to lHinshaw et alj ( l2003h : notice that we 
only include the diagonal terms of the covariance matrix. The final 
estimate of the unresolved sources contribution is the average of 
the spectra from the individual realizations. 

In figure |7] we compare our estimates with WMAP results. 
While our predictions are in good agreement with WMAP5 yr esti- 
mates for the QQ spectra, they show a steeper frequency behaviour. 
In particular, for the final combined power spectrum, we predict 
a residual point source contribution ~ 50% lower than WMAP 
estimates for the source detection strategy discussed above, and 
~ 25% lower if we only mask sources with 5*41 > IJy. Uncer- 
tainty on simulations results is ~ 10%, comparable to the error on 
^0. Notice that in figure |7] we do not plot our estimates of residual 
sources contamination at ^ > 700 for the Q band, as the uncertainty 
in the beam asy mmetries f or that channel introduce significant arti- 
facts at ^ > 500 dHill et a l. 2009). Since we do not use the Q band 
data to estimate the point sources contribution to the final spectra, 
these artifacts do not alter the following discussion. 

The differences between our prediction and WMAP5 esti- 
mates arise from a combination of factors. As discussed above, the 
reference catalog we use to extrapolate high flux sources from 20 
GHz to 94 GHz was given by the subsample of WMAP5 sources 
with with fluxes S23 > 1 Jy, corresponding to the approximate 
completeness limit of WMAP 5yr catalog. WMAP sources with 
5*23 > 1 Jy have an average spectral index a — —0.21 ± 0.30, 
while for sources with 523 < 1 Jy a = —0.09 ± 0.66. The spectral 
index used by the WMAP team in order to derive the amplitude is 
the mean of a Gaussian fit to the spectral index distribution of all 
their detected sources, and is therefore affected by the shallower 
spectral indexes of sources with S23 < 1 Jy Note that we did not 
use the latter set of sources for the simulations, as f or this range 
of flux es and frequencies we adopted the results of I Sadler et al.l 
1 I2OO8 I). These results suggest that sources with fluxes in the 0.2 — 1 
Jy range have an average spectral index a = —0.43 ± 0.31, while 



for the subset of 8 sources with S20 > 0.8Jy,a = — 0.26±0.16, in 
good agreement with estimates for WMAP sources with 5*23 > 1 
Jy. Our simulations reflect this finding. As more information on the 
spectral behaviour of sources with K band flux in the 0.5 — 1.5 
Jy will become available, it will be possible to better establish if 
the low flux sources of WMAP catalog are biasing the estimates of 
the spectral index used for estimating the unresolved point sources 
contribution. 

In addition, WMAP fix the amplitude of the unresolved source 
contribution Aq by a simultaneous fit to the ensemble of cross spec- 
tra obtained using Q, V and W bands, assuming a constant fre- 
quency scaling across the corresponding range of frequencies, 41 
- 94 GHz. Therefore, although the final CMB power spectra are a 
combination of only V and W measurements, the estimated unre- 
solved sources contribution to such spectra depends also on Q data. 
In this work, instead, we estimate the point source contamination 
at a given frequency directly from the point sources maps at that 
frequency, without including information from other bands. In par- 
ticular, our prediction for the final combined VV-i-VW- hWW spec- 
trum d oes not depend on 41 GHz data. As discussed in lNolta et al.l 
( l2009l) . the WMAP5 estimated value of Aq using only V and W 
data is ~ 10% — 40% (depending on the choice of Galactic mask) 
lower than when using Q,V and W data, although the uncertainty in 
the former case is ~ 3 times the uncertainty in the latter. Therefore, 
using information from Q band may lead to a slight overestimate of 
Ao at 60 - 90 GHz. 

As pointed out bv iHuffenberger et al.l ilQOd l2008h the point 
sources correction and the way it is treated in the likelihood func- 
tion affects the determination of the cosmological parameters, in 
particular of the slope of the power spectrum of scalar pertur- 
bations. Us- To test the effect, if any, of our approach on pa- 
rameter es timation, we use COSM OMcQ jLewis & Bridlell2002l) 
-I- PICCEI ( Fendt & Wandeltl |2007|) to run a set of Monte-Carlo 
Markov Chains (MCMC) replacing the source correction term 
in WMAP5 likelihood code with our results. However, we did 
not change the shape of the likelihood function as suggested by 

* http://cosmologist.info/cosmomc/ 
^ http://cosmos.astro.uiuc.edu/pico/ 
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Figure 7. Estimates of the unresolved point sources contribution to CMB power spectra at different WMAP frequencies when applying the masking procedure 
described in the text (solid lines) or simply removing all sources with Q band flux Sai > ^Jy (dashed lines). The dotted lines show a Poisson spectrum with 
normalization Cf^'^ oc Aq with Aq = 0.011 and a = —0.09. The final power WMAP power spectrum is a combination of VV, VW and WW cross 

spectra and does not depend on QQ data. For reference we also show the CMB temperature power spectrum for the best fit WMAP5 cosmological model. 



iHuffenberger et al.l ( l2008h . We consider the six base ACDM pa- 
rameters: physical baryon, ujb, and cold dark matter, ojc, densities; 
optical depth to reionization, Te ; tilt, ris , and amplitude, As of the 
power spectrum of scalar density fluctuations; the angle subtended 
by the sound horizon at recombination, 9. We include marginaliza- 
tion over the amplitude of the Sunayev-Zeldovich effect and lens- 
ing. As expected, since changing the source correction effectively 
alters the shape of the second and third acoustic peaks in the mea- 
sured spectrum, we find that r and, to a lesser degree, 6 and As 
are unaffected. For our default source removal procedure, the shift 
on Us and ujc can reach almost 0.4a, with a noticeable shift also 
for the baryon density, ujt, while if we adopt a more conservative 
masking, the effect can be up to O.Scr. 



7 SUMMARY AND CONCLUSIONS 

We discussed a new method aimed at evaluating the residual 
point source contribution to the WMAP5 power spectrum. We 
implemented a stochastic approach to extrapolate the flux of ra- 
dio sources of the NRAO-VLA Sky Survey (NVSS) to WMAP 
frequencies. Using different point source catalogs with multifre- 
quency information, we divided the 1.4 — 90 GHz in a number of 
frequency steps. At each step an NVSS source is randomly paired 
to a reference source in the catalog covering the relevant frequency 
range, and extrapolated to the upper limit of that range using the 
spectral index of the reference source. The procedure is iterated 



800 times, and the ensemble of simulations provides the probabil- 
ity distribution for the flux each NVSS at the desired frequencies. 

We compared the statistical properties of the simulations with 
those of the sources of the WMAP5 catalog, with other surveys at 
33 GHz and with predictions by previous authors at 94 GHz, find- 
ing an agreement down for fluxes > 0.1 Jy in Ka band, while we 
underestimate the number of sources below ~ 0.1 Jy. The WMAP 
catalog is not complete below S ~ 1 Jy, so point sources with 
S ~ 0.1 Jy account for only a small fraction of the total power 
due to unresolved sources. At higher frequencies predictions are in 
general agreement with previous works. 

An advantage of the method is that it maintains the informa- 
tion about the sources position. In order to test this, we consid- 
ered the sets of point-source-only (i.e. no CMB nor noise) maps 
obtained from NVSS simulations at 41 and 94 GHz, and correlated 
them with both the Q and W band WMAP5 sky maps and with 
mock CMB skies in which the point source population follows the 
same distribution as our model. We found that although the over- 
all level of correlation is low, W band correlations with actual data 
are in good agreement with correlations with the mock data. In Q 
band, instead, the degree of correlation between simulations and ac- 
tual maps is significantly lower than what expected from the mock 
data. In this band foregrounds contamination is significantly higher 
than at higher frequency and may contribute to decreasing the level 
of correlation between simulations and real maps. 

We used the simulated maps to estimate the unresolved source 
correction to WMAP5 measured spectra. As a general result, we 
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Figure 8. Dependence of parameter estimates on the shape of the source 
correction. Curves refer to default WMAP 5yr estimates (black/dotted 
lines), or our prediction when masking sources as described in the text 
(red/soUd) or simply masking sources with > 1 Jy (green/dashed). 



find that unresolved sources in our simulations have a steeper fre- 
quency behaviour than the contribution estimated by the WMAP 
team. This is due to the fact that observation of sources in the 20- 
90 GHz range with S < 1 Jy typically show steeper spectral in- 
dexes than the best fit spectral index from the WMAP5 catalog. 
The observed scaling in the simulations corresponds to an effec- 
tive spectral index a = —0.20 ~ —0.25, compared to WMAP5 
Q = ~ —0.09. Thus, the unresolved point sources contribution 
derived by the simulation either agrees with the WMAP5 estimates 
in Q band but is lower by up to ~ 30 — 40% at 94 GH, depending 
on the way we mask bright sources. The corresponding shift in es- 
timates for parameter like ns and cjc can reach up to ~ 0.3 — OAa. 

Given a sufficient number of simulations, the approach fol- 
lowed in this work allows in principle to estimate the probability 
distribution for the flux at an arbitrary frequency of each source 
from a template low frequency survey. In turn this would allow 
to naturally account for the added uncertainty due to undetected 
sources during the map-making stage, as is done for galactic fore- 
grounds with Gibbs sampling methods. This would require that 
Gibbs samplers can be made to efficiently work at multipoles 
£ ~ 2000 — 3000, where point sources will be relevant for Planck 
and other possibly upcoming high-resolution experiments. Even if 
this were not possible, having an estimate of unresolved source con- 
tributions in each pixel would be helpful in defining a sky mask, e.g. 
by flagging all pixel in which the source contribution is expected to 
exceed a given threshold even if there is no actual detection of a 
point source in that region of the map. 
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